Examine FARM models and their projections v9

workingconservation@gmail.com

Setup libraries

library(tidyverse)
library(terra)
library(MCMCvis)

# ggplot theme set
theme_set(theme_bw())

setup data

wd <- "/Users/elaw/Desktop/LeuphanaProject/BirdModelling/ETH_birds"
results_folder <- paste0(wd, "/Results/")  
version_folder <- "v9/"
mydate <- 20221017
samples <- readRDS(paste0(results_folder, version_folder, "BirdOccMod_dOcc_samples_farm_", mydate, ".RDS"))
data_input <- readRDS(paste0(wd, "/Data/Farm_nimbleOccData_v6_", mydate, ".RDS"))
list2env(data_input[[1]], envir = .GlobalEnv); list2env(data_input[[2]], envir = .GlobalEnv); rm(data_input)

## <environment: R_GlobalEnv>

## <environment: R_GlobalEnv>

rast_folder <- "/Users/elaw/Desktop/LeuphanaProject/ETH_SpatialData/data_10m/Baseline/"
farm_stack <- rast(paste0(rast_folder, "farm_stack_10m.tif"))
pred_stack <- rast(paste0(results_folder, version_folder, "farm_pred_stack.tif"))

rast_vals <- values(pred_stack, na.rm=TRUE) %>% as_tibble()
site_vals <- as_tibble(Xoc); names(site_vals) <- names(rast_vals)

## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
## `.name_repair` is omitted as of tibble 2.0.0.
## ℹ Using compatibility `.name_repair`.

site vars vs raster vars

summarise_all(rast_vals, range) %>% t()

##                [,1]     [,2]
## elevation -2.811507 3.647551
## fldis     -2.126014 3.944279
## sidi1ha   -1.156326 3.160600

summarise_all(site_vals, range) %>% t()

##                [,1]     [,2]
## elevation -1.788999 1.629391
## fldis     -1.380924 2.606278
## sidi1ha   -1.156326 2.282887

Plot histograms

Site sampled variables (green) are typically a biased sub-section of the range of variables in the rasters (purple). This is not as much of an issue as in forest.

plot_hists <- function(myvar, binwidth=0.25){
  mc <- viridis::viridis(3)
  ggplot(rast_vals, aes(x={{myvar}}, y = stat(density))) +
    geom_histogram(binwidth = binwidth, fill=mc[1], alpha = 0.5) + 
    geom_histogram(data = site_vals, binwidth = binwidth, fill=mc[2], alpha = 0.5)
}

plot_hists(elevation)

plot_hists(fldis)

plot_hists(sidi1ha)

Examine betalpsi parameters from each species

# select a reasonable set of samples 
mydraw <- 501:1500
mysamples <- list(
  chain1 = samples$chain1[mydraw, ],
  chain2 = samples$chain2[mydraw, ],
  chain3 = samples$chain3[mydraw, ]
)

# extract betas
b0 <- MCMCsummary(mysamples, 
                  params = 'lpsi\\[\\d{1,3}\\]', ISB = FALSE, exact = FALSE,  
                  round = 2)
b1 <- MCMCsummary(mysamples, 
                  params = 'betalpsi\\[\\d{1,3}[,][ ][1]\\]', ISB = FALSE, exact = FALSE,  
                  round = 2)
b2 <- MCMCsummary(mysamples, 
                  params = 'betalpsi\\[\\d{1,3}[,][ ][2]\\]', ISB = FALSE, exact = FALSE,  
                  round = 2)
b3 <- MCMCsummary(mysamples, 
                  params = 'betalpsi\\[\\d{1,3}[,][ ][3]\\]', ISB = FALSE, exact = FALSE,  
                  round = 2)

tibble(
  ElevationMean = summary(b1$mean > 0), ElevationMedian = summary(b1$`50%`>0),
  FarmDistanceMean = summary(b2$mean > 0), FarmDistanceMedian = summary(b2$`50%`>0),
  Sidi1haMean = summary(b3$mean > 0), Sidi1haMedian = summary(b3$`50%`>0)
) %>% t %>% 
  as_tibble(rownames = "SummaryType") %>% 
  select(-V1, Negative=V2, Positive=V3)

## # A tibble: 6 × 3
##   SummaryType        Negative Positive
##   <chr>              <chr>    <chr>   
## 1 ElevationMean      135      41      
## 2 ElevationMedian    138      38      
## 3 FarmDistanceMean   157      19      
## 4 FarmDistanceMedian 157      19      
## 5 Sidi1haMean        22       154     
## 6 Sidi1haMedian      22       154

Plot the betalpsi params

plotbeta <- function(bp, bname, fgroup=NULL, fgroupName){
  bp <- bind_cols(bp, 
                  fspp=c(fspp, rep("unknown", Nadd)), 
                  mig=c(mig, rep("unknown", Nadd)),
                  fnDiet=c(fnDiet, rep("unknown", Nadd)),
                  invDiet=c(invDiet, rep("unknown", Nadd)))
  ggplot(arrange(bp,`50%`), 
         aes(x=`50%`, xmin=`2.5%`, xmax=`97.5%`,
             y=1:M))+
    geom_linerange(aes(color={{fgroup}}))+
    geom_point(aes(x=mean), color="darkgrey")+
    geom_point(aes(color={{fgroup}}))+
    geom_vline(xintercept=0, color="red")+
    scale_color_viridis_d()+
    labs(x=paste0(bname," beta parameter:\n95% HDI (lines), mean (grey points), median (color points)"),
         y="Species, sorted by median beta",
         color = fgroupName)+
    theme(axis.text.y=element_blank(),
          axis.ticks.y=element_blank())
}

plotbeta(b1, "Elevation", fspp, "Forest species")

plotbeta(b2, "Farm distance", fspp, "Forest species")

plotbeta(b3, "Landscape diversity (sidi1ha)", fspp, "Forest species")

plotbeta(b1, "Elevation", mig, "Migratory species")

plotbeta(b2, "Farm distance", mig, "Migratory species")

plotbeta(b3, "Landscape diversity (sidi1ha)", mig, "Migratory species")

Plot beta relationships

newdata <- tibble(
  elevation = seq(-2.811507, 3.647551, length.out=100) %>% rep(2),
  mElevation = mean(rast_vals$elevation),
  fldis = seq(-2.126014, 3.944279, length.out=100) %>% rep(2),
  mFldis = mean(rast_vals$fldis),
  sidi1ha = seq(-1.156326, 3.160600, length.out=100) %>% rep(2),
  mSidi1ha = mean(rast_vals$sidi1ha)
)

elev_mfd_mld_psi <- matrix(data=NA, nrow = nrow(newdata), ncol = M)
for(k in 1:M){
  logitpsi <- b0$mean[k]  +
    b1$mean[k] * newdata$elevation + 
    b2$mean[k] * newdata$mFldis +
    b3$mean[k] * newdata$mSidi1ha 
  
  elev_mfd_mld_psi[,k] <- exp(logitpsi)/(exp(logitpsi) + 1)
}

melev_fldis_mld_psi <- matrix(data=NA, nrow = nrow(newdata), ncol = M)
for(k in 1:M){
  logitpsi <- b0$mean[k]  +
    b1$mean[k] * newdata$mElevation + 
    b2$mean[k] * newdata$fldis +
    b3$mean[k] * newdata$mSidi1ha 
  
  melev_fldis_mld_psi[,k] <- exp(logitpsi)/(exp(logitpsi) + 1)
}

melev_mfd_sidi_psi <- matrix(data=NA, nrow = nrow(newdata), ncol = M)
for(k in 1:M){
  logitpsi <- b0$mean[k]  +
    b1$mean[k] * newdata$mElevation + 
    b2$mean[k] * newdata$mFldis +
    b3$mean[k] * newdata$sidi1ha 
  
  melev_mfd_sidi_psi[,k] <- exp(logitpsi)/(exp(logitpsi) + 1)
}


out_psi <- elev_mfd_mld_psi %>% 
  as_tibble() %>% 
  bind_cols(newdata) %>% 
  pivot_longer(cols = V1:V176, names_to = "species", values_to = "elev_mfd_mld_psi") %>% 
  mutate(species = str_replace(species, "V", "sp")) %>% 
  bind_cols(
    melev_fldis_mld_psi %>% 
      as_tibble() %>% 
      pivot_longer(cols = V1:V176, names_to = "species", values_to = "melev_fldis_mld_psi") %>% 
      select(-species)
  ) %>% 
    bind_cols(
    melev_mfd_sidi_psi %>% 
      as_tibble() %>% 
      pivot_longer(cols = V1:V176, names_to = "species", values_to = "melev_mfd_sidi_psi") %>% 
      select(-species)
  ) %>% 
  add_column(
    b1 = rep(b1$mean, 200),
    b2 = rep(b2$mean, 200),
    b3 = rep(b3$mean, 200),
    fspp = rep(c(fspp, rep(FALSE, Nadd)),200), 
    fnDiet=rep(c(fnDiet, rep(FALSE, Nadd)),200), 
    invDiet=rep(c(invDiet, rep(FALSE, Nadd)),200),
    observed=rep(c(rep("observed", Nobs), rep("absent", Nadd)), 200)
  ) %>% 
  mutate(fspp = if_else(fspp==TRUE, "fspp", "other"),
         fnDiet = if_else(fnDiet==TRUE, "fnDiet", "other"),
         invDiet = if_else(invDiet==TRUE, "invDiet", "other"))

plotpsi <- function(xvar, yvar, xname, yname, oxmin=-Inf, oxmax=Inf){
  ggplot(out_psi, 
         aes(x={{xvar}}, y={{yvar}}, color=species)) +
    # add rectangles showing non-sampled projection area
    annotate("rect", xmin = -Inf, xmax = oxmin, ymin = -Inf, ymax = Inf, 
             alpha =0.3)+
    annotate("rect", xmin = oxmax, xmax = Inf, ymin = -Inf, ymax = Inf, 
             alpha =0.3)+
    geom_line(alpha=0.5) + 
    scale_color_viridis_d()+
    theme(legend.position = "none") +
    labs(x=xname, y=yname)
}

plotpsi(elevation, elev_mfd_mld_psi, 
        "Elevation (center scaled, Farm distance and sidi1ha at mean)", "psi", 
        -1.788999, 1.629391)

plotpsi(fldis, melev_fldis_mld_psi, 
        "Farm distance (center scaled, Elevation and sidi1ha at mean)", "psi", 
        -1.380924, 2.606278)

plotpsi(sidi1ha, melev_mfd_sidi_psi, 
        "Landscape diversity (sidi1ha center scaled, Elevation and farm distance at mean)", "psi", 
        -1.156326, 2.282887)

Sum species over the parameters

plotSR <- function(xvar, y, xname, yname, oxmin=-Inf, oxmax=Inf){
  bind_cols(
    newdata,
    SR = rowSums(y),
    SRfspp = rowSums(y[,fspp]),
    SRmig = rowSums(y[,mig]),
    SRfnDiet = rowSums(y[,fnDiet]),
    SRinvDiet = rowSums(y[,invDiet])
  ) %>% 
    pivot_longer(cols = SR:SRinvDiet, names_to = "SpeciesSelection", values_to = "SR") %>% 
  ggplot(., 
         aes(x={{xvar}}, y=SR, color=SpeciesSelection), lwd=2) +
    # add rectangles showing non-sampled projection area
    annotate("rect", xmin = -Inf, xmax = oxmin, ymin = -Inf, ymax = Inf, 
             alpha =0.3)+
    annotate("rect", xmin = oxmax, xmax = Inf, ymin = -Inf, ymax = Inf, 
             alpha =0.3)+
    geom_line() + 
    scale_color_viridis_d()+
    labs(x=xname, y=yname)
}

plotSR(elevation, elev_mfd_mld_psi, 
        "Elevation (center scaled, Farm distance and sidi1ha at mean)", "Species Richness (sum PSI)", 
        -1.788999, 1.629391)

plotSR(fldis, melev_fldis_mld_psi, 
        "Farm distance (center scaled, Elevation and sidi1ha at mean)", "Species Richness (sum PSI)", 
        -1.380924, 2.606278)

plotSR(sidi1ha, melev_mfd_sidi_psi, 
        "Landscape diversity (sidi1ha center scaled, Elevation and farm distance at mean)", "Species Richness (sum PSI)", 
        -1.156326, 2.282887)

Probability of detection

MCMCsummary(mysamples, 
                  params = "betalp",  
                  round = 2)

##            mean   sd  2.5%   50% 97.5% Rhat n.eff
## betalp[1]  0.11 0.04  0.04  0.11  0.18    1  2458
## betalp[2]  0.05 0.04 -0.01  0.05  0.12    1  3066
## betalp[3] -0.20 0.04 -0.28 -0.20 -0.13    1  2929

# extract betas for detection
d0 <- MCMCsummary(mysamples, 
                  params = "lp",  
                  round = 2)
d1 <- MCMCsummary(mysamples, 
                  params = 'betalp\\[[1]\\]', ISB = FALSE, exact = FALSE, 
                  round = 2)
d2 <- MCMCsummary(mysamples, 
                  params = 'betalp\\[[2]\\]', ISB = FALSE, exact = FALSE,  
                  round = 2)
d3 <- MCMCsummary(mysamples, 
                  params = 'betalp\\[[3]\\]', ISB = FALSE, exact = FALSE,  
                  round = 2)

logitp_v1 <- matrix(d0$mean, length(d0$mean), dim(Xob)[1])  +
  d1$mean * Xob[,1,1] + 
  d2$mean * Xob[,1,2] +
  d3$mean * Xob[,1,3] 
logitp_v2 <- matrix(d0$mean, length(d0$mean), dim(Xob)[1])  +
  d1$mean * Xob[,2,1] + 
  d2$mean * Xob[,2,2] +
  d3$mean * Xob[,2,3] 

p_v1 <- exp(logitp_v1)/(exp(logitp_v1) + 1)
p_v2 <- exp(logitp_v2)/(exp(logitp_v2) + 1) 

# mean probability of detection, all species, all sites, on visit 1 and visit 2
mean(p_v1); mean(p_v2)

## [1] 0.05202893

## [1] 0.06399845

# plot probability of detection
tibble(
  mean_p_v1 = rowMeans(p_v1),
  mean_p_v2 = rowMeans(p_v2)
) %>% 
  arrange(-mean_p_v2, -mean_p_v1) %>% 
  add_column(species = 1:M) %>% 
  ggplot(., aes(x=species, ymin=mean_p_v1, ymax=mean_p_v2)) +
  geom_hline(yintercept=mean(p_v1), color="red", lty="dotted") +
  geom_hline(yintercept=mean(p_v2), color="red") +
  geom_linerange() +
  geom_point(aes(y=mean_p_v1), pch=1) +
  geom_point(aes(y=mean_p_v2), color="black") +
  labs(x="Species, sorted by probability of detection (visit 2)", 
       y="Mean probability of detection over all sites\nOpen circles = visit 1; Closed circles = visit 2")